Simulations of Disordered Bosons on Hyper-Cubic Lattices 
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Abstract 

We address computational issues relevant to the study of 
disordered quantum mechanical systems at very low tem- 
peratures. As an example we consider the disordered Bose- 
Hubbard model in three dimensions directly at the Bose- 
glass to superfluid phase transition. The universal aspects 
of the critical behaviour are captured by a (3 + 1) dimen- 
sional link-current model for which an efficient 'worm' al- 
gorithm is known. We present a calculation of the distribu- 
tion of the superfluid stiffness over the disorder realizations, 
outline a number of important considerations for perform- 
ing such estimates, and suggest a modification of the link- 
current Hamiltonian that improves the numerical efficiency 
of the averaging procedure without changing the universal 
properties of the model. 



1. Introduction 

The Bose-Hubbard model was first studied in the con- 
text of liquid helium in a disordered medium |6|. Interest 
in the model has recently grown with the progress achieved 
in trapping ultra-cold atomic gases in optical lattice poten- 
tials. The model describes the competition between tun- 
neling and on-site interactions in a lattice of bosons. It dis- 
plays several zero-temperature quantum phases that are now 
clearly attainable in the laboratory |7|. Notably, a localized 
Mott insulating (MI) phase exists when the tunneling be- 
tween sites is small, while at higher tunneling the system 
becomes a superfluid (SF). In the presence of disorder an- 
other localized phase, the Bose-glass (BG), exists between 
the Mott insulator and the superfluid. In the present study 
we focus on the Bose-glass to superfluid transition since it 
exists only at finite disorder and thus provides a clear ex- 
ample of a quantum phase transition for which disorder is 
relevant. Although this model has been extensively studied 
in one- and two-dimensions, the nature of the phase transi- 
tion in three and higher dimensions has received relatively 
little attention. 

Scaling theories based on generalized Josephson rela- 



tions and the finite compressibility of the superfluid and 
Bose-glass phases indicate that the dynamic correlation ex- 
ponent z is equal to the number of spatial dimensions 1 6 1 . 
This feature, which is supported by analytical and numeri- 
cal arguments 1 8 9 , 13 1 in low dimensions suggests that the 
model has an unusual approach to mean-field behaviour — 
and may not have an upper critical dimension at all — 
invalidating standard renormalization group approaches. 
Numerical work above two dimensions is difficult due to 
the large volumes of the systems, and the algorithmic slow 
down of the Monte Carlo averaging procedure. 

Of particular interest are the distributions of thermody- 
namic observables over the disorder. They are typically far 
from Gaussian in nature, rendering standard estimates of 
statistical error invalid for smaller sample sizes and necessi- 
tating calculations for a large number of disorder configura- 
tions. Moreover, the behaviour of these distributions for in- 
creasing lattice size is directly relevant to the break-down of 
self-averaging 1 2 1, the quantum Harris Criterion 1 5 1 1 1, and 
the effect of disorder on quantum critical phenomena 1 12 1. 
More efficient ways of performing disorder averages has 
also been proposed 1 3 1. These latter developments are how- 
ever too computationally demanding for the present model. 

The universal properties of the rf-dimensional Bose- 
Hubbard model are captured by a (d + l)-dimensional 
classical link-current representation for which an efficient 
worm-like Monte Carlo algorithm exists. While the worm 
algorithm represents a drastic improvement over earlier, 
Metropolis-like algorithms, the computational demands in- 
crease dramatically in higher dimensions, limiting the pre- 
cision of the numerical analysis. Since the system is dis- 
ordered, the calculations involve performing many Monte 
Carlo simulations of the system at the same parameters with 
different realizations of the disorder. Hence, the calcula- 
tions are very well suited for parallelization and a linear 
speed up can be achieved with a relatively modest program- 
ming effort. Without such a linear speed up the calcula- 
tions we report on would have been almost impossible. Par- 
allelization is performed straightforwardly with MPI: each 
processor performs Monte Carlo simulations serially for a 
given disorder realization, then the results are collected and 
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written to disk. As many as several thousand disorder re- 
alizations need to be performed at each point in parameter 
space, and each individual Monte Carlo simulation can take 
up to three hours, depending on the system size and the pa- 
rameters of the model. The length of time required to per- 
form one simulation for a given disorder realization dictates 
how large a system we can reasonably study and so it is es- 
sential to consider carefully how much computational effort 
to invest in each such simulation. 

Numerical study of this transition in three dimensions 
presented a number of difficulties in calculating the disor- 
der distributions and their averages. The main results of the 
study will be presented elsewhere [ 10 1, but we outline here 
the procedure that was used to estimate these distributions, 
and suggest a modification to the link-current Hamiltonian 
which improves the efficiency of these numerical estimates. 
The remainder of the introduction discusses the link-current 
Hamiltonian and the finite-size scaling theory on which our 
numerical approach relies. Section|2]describes how we en- 
sure that the simulation of each disorder realization has been 
properly equilibrated. Section [3] describes a modification 
of the link-current model that improves the numerical ef- 
ficiency of each simulation without affecting the universal 
details. We then conclude with some remarks about the gen- 
eral applicability of these considerations. 

1.1. Model and Scaling Theory 

The Bose-Hubbard Hamiltonian, including an on-site 
disorder in the chemical potential is 1 6 1 
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The second quantized boson operators describe a tunneling 
process coupled by t and an on-site, repulsive interaction 
U, on a hyper-cubic lattice. The disordered chemical po- 
tential /i r is distributed uniformly on [fj, — A, [i + A] so that 
A controls the strength of the disorder. At finite disorder, 
the system undergoes a phase transition from a Bose-glass 
insulating phase (low t, high U) to a superfiuid phase (high 
t, low U). The model can be transformed via the quantum 
rotor model to the (d + 1) link-current model 1 13 1. The link- 
current Hamiltonian is given by 
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The integer currents J( r T ) are situated on the bonds of the 
lattice and obey a divergenceless constraint 



The resulting loops are interpreted as currents of bosons 
hopping about on the lattice (specifically they are fluctua- 
tions from an average density, so that the currents are per- 
mitted to be negative). The coupling K controls the ratio 
between t and U : at low K the interaction U dominates 
and the system is insulating, while at high K the tunneling 
dominates and the system condenses into a superfiuid. 

The two phases can be distinguished by the superfiuid 
stiffness, p, which is proportional to the superfiuid density. 
The stiffness is defined as the response of the free energy to 
a twist in the boundary conditions and is indicative of global 
phase coherence. In the link-current model, the stiffness 
is proportional to the square of the winding number in the 
spatial dimensions: 

P= JJ^T t [<»»>] W 

The angle brackets (•) denote a thermal average, while the 
square brackets [•]„„ denote an average over disorder real- 
izations. The winding numbers (n 7 = L~ 1 J2 r r ^rr f° r 
7 = x,y, z, t) of the lattice in each direction are just the 
number of current loops that have wound all the way about 
the periodic lattice. They are always integers. 

Dynamics and statics are both essential to the critical 
behaviour of quantum phase transitions. They are char- 
acterized by independent spatial and temporal correlation 
lengths (£ and £ r ) which define the dynamic critical expo- 
nent z: 

(*-")*. 8 — ~~~~r7~~~~ ■ (5) 

Here v is the correlation length exponent. This implies 
that quantities at the critical point scale as a function of 
two arguments. The superfiuid stiffness (which diverges as 
p ~ £ d + z ~ 2 at K c ) thus scales as 

1 
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If we hold the second argument constant, the critical point 
can be located by plotting curves of pL d+z ~ 2 for various 
linear system sizes L. Since 6 = at K c , K c will be the 
value of K at which these curves intersect. This unfortu- 
nately requires that we guess at the value of z before we 
begin. In principle, we are free to set the aspect ratio 



L T /L Z 



(7) 



as we see fit to find the critical point; in practice however, as 
we discuss below, the numerics work better near an optimal 
aspect ratio where (£,/L) z ~ £, T /L T implying a — 0(1). 

2. Equilibration 



f—x 7 y : z : T 



J (r,r) 



0. 



(3) When considering disordered systems, the average 

[(■)}av of an arbitrary observable (denoted by •) such as the 
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stiffness p must be calculated over a whole set of disorder 
configurations, performing independent Monte Carlo simu- 
lations on each particular realization. We must then decide 
how many Monte Carlo sweeps (i„) to perform on each sim- 
ulation, and how many disorder realizations (A^) to aver- 
age these over. There are correspondingly two sources of 
statistical error [4 1: the error 5tP in the estimate of the ther- 
mal average: 

{p)=p + 5 T p, (8) 
and the overall error <5[p] au in the disorder average: 



[{P)]av = [P + 5 T p]aV = [p]av + S\p] c 



(9) 



It is of particular interest to review [4 1 how to correctly ob- 
tains the disorder average of the square of a thermodynamic 
observable such as the energy. Such a quantity would be 
needed for calculating for instance the specific heat, Cy. In 
this case we have for a single disorder realization: 

{E)=E + S T E. (10) 

If we now want to calculate the disorder average [(E) 2 ] av 
we encounter a slight problem: 



[(StE) 2 



2E[S T E] av . (11) 



It is natural to assume that [8TE] av will yield zero when a 
sufficiently large number of disorder realizations are used. 
However, [(8TE) 2 ] av will be non zero and will yield a sys- 
tematic error unless infinitely precise thermal averages can 
be obtained for each disorder realization. In order to cir- 
cumvent this problem and correctly calculate such a dis- 
order average, one can run 2 independent simulations, re- 
ferred to as "replicas", of a given disorder realization |4|. 
We denote them by a and (3. The above disorder average 
should then be calculated in the following way: 



[(E)\ v = [(E a )(E^)] av 

= [E 2 ] av + [8 T E a 5 T E% v 
+E[5 T E a ] av + E[8 T E^\ 



(12) 



The term [6 T E a S T E% v will now also correctly average to 
zero since the thermal errors from each replica are indepen- 
dent random variables. In our calculations we always run at 
least two independent replicas of a given disorder realiza- 
tion with the goal of correctly calculating averages as out- 
lined above. For higher powers of thermal averages more 
replicas are needed. As we shall see below, having several 
independent runs of a given disorder realization allow also 
for very convenient and indispensable checks of the equili- 
bration of the calculations. 

It is important to note that even for very large system 
sizes, the average should be taken over as many disorder 
realizations as possible. One might assume that for large 
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Figure 1. Hamming distances on an 
8x8x8x64 lattice at K c = 0.19 calculated over 
a set of 1000 disorder realizations with 
t = 3 x 10 7 . For the Hamming distances be- 
tween the two replicas a and (3, t is the total 
number of Monte Carlo sweeps performed. 
For the Hamming distances between replica 
a and its configuration a at t , t is the num- 
ber of sweeps performed after the initial t 
sweeps have been performed. The conver- 
gence of the curves indicates t r « 3 x io 5 . 



system sizes a smaller set of disorder realizations would 
suffice, since the properties of the system will average out 
spatially. For many disordered systems this assumption is 
false — self-averaging breaks down 1 1 1. Even in the thermo- 
dynamic limit (where in principle one could find any par- 
ticular finite disorder realization somewhere in the infinite 
system) the width of the distribution of P(p) remains finite. 
Moreover, such distributions are typically far from Gaus- 
sian (they often have a particularly 'fat' tail). Error esti- 
mates based on Gaussian distributions are thus only valid 
for very large No- The best approach is then to spend a 
minimal amount of computational time on each realization, 
and then rely on the disorder average to control the statisti- 
cal errors. As usual, however, each Monte Carlo simulation 
must be properly equilibrated to ensure that (p) is unbiased. 
Since the lattice starts in an artificial (and non-equilibrium) 
configuration, we throw out the first to of the t n sweeps be- 
fore sampling the remaining t s = t n — to configurations. As 
long as to is greater than the relaxation time of the algorithm 
t r , we sample only equilibrium configurations, and our es- 
timates of the thermal averages will be unbiased. Since we 
do not want to spend all of our computational efforts equi- 
librating systems, each disorder realization is typically run 
for an additional to steps once it has reached equilibrium. 
To confirm that we have chosen to greater than t r , we 



3 



perform each simulation independently on two replicas with 
different initial configurations. We can then define 'Ham- 
ming distances' 1 13 1 between the two replicas a and (3 after 
performing t Monte Carlo sweeps on their initial configu- 
ration, and between the replica a at sweep t + to and its 
configuration ao at sweep to when sampling begins: 

T (r,r) 

Kj- T (t) = j)r^ [ J -(^)W-4(r,.)W] 2 - (") 

T (r,r) 

These measures are then averaged over the disorder realiza- 
tions. At the beginning of each simulation, the initial con- 
figurations of a and (3 are quite different. Hence H a ^(t) 
will be large for small t, but diminish as the two configu- 
rations equilibrate. On the other hand, shortly after to, the 
configuration a will not have changed substantially from 
ao, so Ha, >OCo {t) will be initially small, but increase as more 
sweeps are performed. If to has in fact been chosen greater 
than t r , these two measures will converge on the same value 
in t r sweeps at which point the configurations of a, ao, 
and (3 will be independently drawn from the same popula- 
tion of equilibrium configurations. Figure[2shows the two 
Hamming distances plotted as a function of t at K c , aver- 
aged over Nd = 1000 disorder realizations. Results are 
shown for Hamming distances defined in terms of spatial 
(x, y, z) currents (open symbols) and the temporal (r) cur- 
rent (filled symbols). Both the spatial and temporal Ham- 
ming distances converge, indicating that t r ks 3 x 10 5 for 
an 8x8x8x64 lattice at the critical point. 

From these observations it would seem reasonable to 
generate t s = 3 x 10 5 configurations after equilibration at 
each disorder realization in order to calculate the disorder 
average. However, if we look more closely at the distribu- 
tion P((p)) generated by using t s = 3 x 10 5 configurations, 
we find a large peak in the distribution at (p) = (see the 
uppermost graph of Fig. 0. For many disorder realizations, 
the Monte Carlo algorithm never generates an equilibrium 
configuration with a non-zero winding number (this can be 
verified by looking directly at the data set). If we generate 
more configurations at equilibrium for each disorder real- 
ization (that is, we increase t s ), the shape of the distribu- 
tion changes — the peak at (p) = shrinks and a broader 
one grows at a non-zero value. After one to three million 
sweeps, the peak at (p) = disappears and the distribution 
stops evolving. This effect is a result of the discrete na- 
ture of the winding number, n x , and the strong correlation 
between configurations generated by a successive Monte 
Carlo sweeps. Since the superfluid stiffness is defined in 
terms of n x , it is difficult to numerically resolve the dif- 
ference between a superfluid stiffness of zero and a small 
fractional value. For instance, a disorder realization with 
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Figure 2. Evolution of the distribution of 
(pi 4 ) as a function of t s (the number of sam- 
ples gathered at equilibrium for each disor- 
der realization) for the same sample disorder 
realizations used in Fig . 03 The vertical axes 
are offset for clarity. The peak at (pi 4 ) = o 
persists for t s > t r , giving way to a broader 
peak at a non-zero (pL A ). The distribution 
continues to change until t s «3x 10 6 upon 
which the distribution stops evolving. 
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Figure 3. Distribution of (pL 4 ) calculated 
with different t but the same t s over two dif- 
ferent sets of 1000 disorder realizations. The 
distribution is unchanged. 



pL 4 — 0.05 corresponds to an average squared winding 
number of (n|) = jj(pL 4 ) — 1/160, roughly implying 
that only one out of every 160 independent configurations 
has a non-zero winding number. Since the algorithm gener- 
ates approximately one fully independent configuration ev- 
ery t r sweeps, at t s — 3 x 10 5 one would expect most runs to 
find (pL 4 ) = 0. (Here we assume that t r is proportional to 
the autocorrelation time r.) However, if the average is taken 
over a further 10 2 independent configurations, the estimate 
will converge on the small but finite true average. Resolving 
the true shape of the distribution thus requires much longer 
runs at each realization of the disorder. That the lattice is 
in fact equilibrated at t r as determined by the convergence 
of the Hamming distances is supported by the fact one can 
generate the same distributions independent of to (so long 
as to > t r ), Figure[3]shows the distribution of (pL 4 ) calcu- 
lated over two different sets of 1000 disorder realizations, 
one with i = 3x 10 5 , the other with t = 3 x 10 7 . In 
both cases t s = 3 x 10 7 configurations were generated after 
equilibration for each disorder realization. 



An immediate conclusion to be drawn from these results 
is that, if attention is not paid to either the Hamming dis- 
tances or the convergence of the complete distribution of 
the thermodynamic observables over the disorder, then one 
is very likely to be mislead by the results obtained. In par- 
ticular, it is clear from the results in Figure |2] that if t s or 
Nd are too small, then [{■)]„„ will be too small and the er- 
ror bars will also be misleadingly small. As t s and No are 
increased, [(•)]„„ will also likely increase due to the tails of 
the distribution while at the same time the associated error 
bars might very well remain roughly constant. 



3. Coupling Anisotropy 

In the context of the Bose-glass to superfluid transition, 
one approach to improving the efficiency of the calculation 
is to increase the numerical value of the stiffness at the crit- 
ical point. One possible means of achieving this is to in- 
crease the aspect ratio a between the spatial and temporal 
lattice sizes. In a finite system, the correlation length can be 
no longer than the size of the system; but since the temporal 
and spatial correlation length are related, a small a that trun- 
cates £ T will in turn restrict £ to be smaller than the spatial 
extent of the lattice. Since a global current loop depends 
on correlation across the whole length of the lattice, this 
will limit the number of winding events and in turn reduce 
p. At the critical point, the value of \(pL 4 )] av should in- 
crease with a and it seems natural to assume that it will con- 
tinue to increase at least up until £, T /L T ~ (£/L) z , where 
a = 0(1). Simulations are then best performed at this op- 
timal aspect ratio; unfortunately they prove too large for us 
to simulate. The problem here is the relatively small prob- 
ability of generating configurations with non-zero winding 
number in the spatial direction if a very small a is used. 
One could then ask if equivalent models can be found with 
higher winding numbers at K c for the same aspect ratio. 
Such models would likely have an optimal a at smaller val- 
ues than the present one thereby decreasing the numerical 
overhead. It turns out that for the present model we can 
exploit the universality of the critical behavior to arrive at 
equivalent models that indeed have this desired property. 

We now outline how we arrive at these equivalent mod- 
els. The transformation from the Bose-Hubbard Hamil- 
tonian to the effective link-current Hamiltonian described 
in Ref. 1131 yields, near the end of the procedure, an 
anisotropic action 



J ^ (r,r) 



(14) 



where At is the width of each time slice, J' = 

E„=x,y,z j2 > and 

K t = UAr, K x = -2rn(-tAr/2). (15) 

The simplest approach here is to set KK X = KK T = 
1, which yields the isotropic link-current model as stated 
above with /i r = (E T /U. There is a freedom here, however: 
one can instead set KK X = 1 and KK T = 7, and introduce 
an anisotropy between the space and time couplings in the 
link-current model without affecting the universal details: 
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(16) 
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Figure 4. Crossings of [(pL A )] av for three 
values of the anisotropy 7 = 0.3, 1, and 3. 
The critical point increases as a function of 
7, as does the value of [{pL A )] av at the criti- 
cal point. 



The link-current couplings can be mapped back to the 
Bose-Hubbard tunneling and on-site disorder using dl5l >: 

^ = 2 7 (Ke^ K y\ (17) 

If the transition occurs at a particular ratio U/t, an increase 
in 7 implies then an increase in K c , for K c < 1/2. More 
importantly, increasing 7 freezes out the temporal dynam- 
ics of the link-current model. Both of these effects should 
speed up the spatial dynamics of the 'worm' algorithm, and 
since p is defined in terms of these winding numbers, by 
changing 7 we can tune the value of pL A at the critical point. 
Figure^shows crossings in pL 4 for 7 = 0.3, 1, and 3. As 
expected, increasing 7 increases the value of K c and the 
magnitude of the stiffness at the crossing, hence improving 
the sampling efficiency of the algorithm. Note that this has 
been achieved without changing the aspect ratio, a. 

Due to space constraints we do not show results for the 
critical exponents at different 7. However, we have studied 
them in detail and they are indeed independent of 7 as one 
would expect from universality arguments. We note that 
it would be of considerable interest to study the evolution 
of the distribution of thermodynamic variable, as shown in 
Fig.|2]for 7 = 1, for different values of 7. Preliminary re- 
sults indicates that the width of the distributions vary with 
7, likely increasing monotonically with this parameter. Due 



to time constraints we leave a detailed investigation for fu- 
ture study. 

4. Conclusions 

Numerical studies of disordered systems are notoriously 
difficult due to the large amount of computational resources 
required and to the many subtle sources of systematic error. 
The procedure we have presented to estimate disorder distri- 
butions and their averages highlights some of the potential 
pitfalls. Under-sampling the disorder distribution (setting 
too low), under-sampling each individual distribution 
(setting t s too low), and failing to properly equilibrate each 
simulation (setting to too low) can all lead to erroneous es- 
timates. The two procedures presented above provide con- 
firmation that these pitfalls have been avoided. The conver- 
gence of the Hamming distances provides a measure of t r 
and confirms that the Monte Carlo simulations are generat- 
ing equilibrium configurations. The disorder distributions 
themselves can contain artifacts of the Monte Carlo averag- 
ing procedure; they should not demonstrate any dependence 
on to or t s if these have been set sufficiently large. These 
procedures can be generalized to other systems. 

In the context of the Bose-Hubbard model, the discrete 
nature of the winding number in the link-current represen- 
tation forces the investment of a large amount of computa- 
tional resources in simulating each individual disorder re- 
alization (t s must be set much greater than t r ). However, 
by adjusting the anisotropy between the spatial and tem- 
poral coupling strengths in the link-current model, we can 
increase the magnitude of the superfluid stiffness at the crit- 
ical point. This improves the numerical efficiency of the 
calculation without affecting the universal details of the crit- 
ical behaviour. This observation may prove useful for future 
Monte Carlo studies of the Bose-Hubbard model in the link- 
current representation. 
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